MODULE fefmm_LIP6_Fij

   USE fmm_global_paras
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: fefmm_build_Fij_LIP6

   REAL(REALK),   PARAMETER :: T_VALUE_THR = 0.9D0
   INTEGER(INTK), PARAMETER :: LIPN = 6

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_build_Fij_LIP6(FEdim,NFS,FEgrain,t,Fij)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: FEdim, NFS
      REAL(REALK),   INTENT(IN)  :: FEgrain, t
      REAL(REALK),   INTENT(OUT) :: Fij(:,:)

      IF ( t > T_VALUE_THR ) THEN
         CALL fefmm_build_Fij_LIP6_analytic(FEdim,NFS,FEgrain,t,Fij)
      ELSE
         CALL fefmm_build_Fij_LIP6_numerical(FEdim,NFS,FEgrain,t,Fij)
      END IF

   END SUBROUTINE fefmm_build_Fij_LIP6

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_build_Fij_LIP6_analytic(FEdim,NFS,FEgrain,t,Fij)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: FEdim, NFS
      REAL(REALK),   INTENT(IN)  :: FEgrain, t
      REAL(REALK),   INTENT(OUT) :: Fij(:,:)

      INTEGER(INTK) :: i,k, u,v

      Fij = zero

      DO i = 1, FEdim 
         v = i + FEdim*((NFS-1)/2)  ! V-point

         DO k = 1, NFS*( (FEdim -1)/LIPN )   ! finite elements
            u = 1 + LIPN*(k-1)    ! basic rho-point

            ! 1st LIP
            Fij(u,i)   = Fij(u,i)   + get_F1_LIP6(t,FEgrain,v,u)
            ! 2nd LIP
            Fij(u+1,i) = Fij(u+1,i) + get_F2_LIP6(t,FEgrain,v,u)
            ! 3rd LIP
            Fij(u+2,i) = Fij(u+2,i) + get_F3_LIP6(t,FEgrain,v,u)
            ! 4th LIP
            Fij(u+3,i) = Fij(u+3,i) + get_F4_LIP6(t,FEgrain,v,u)
            ! 5th LIP
            Fij(u+4,i) = Fij(u+4,i) + get_F5_LIP6(t,FEgrain,v,u)
            ! 6th LIP
            Fij(u+5,i) = Fij(u+5,i) + get_F6_LIP6(t,FEgrain,v,u)
            ! 7th LIP
            Fij(u+6,i) = Fij(u+6,i) + get_F7_LIP6(t,FEgrain,v,u)
!if (abs(Fij(u,i)) > 1d-5) print *, "ana_Fij: ", u,i, v-u, get_F1_LIP6(t,FEgrain,v,u)

         END DO

      END DO

   END SUBROUTINE fefmm_build_Fij_LIP6_analytic

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_build_Fij_LIP6_numerical(FEdim,NFS,lamb,t,Fij)

      USE fefmm_quadrature, ONLY: fefmm_get_roots_and_weights

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN)  :: FEdim, NFS
      REAL(REALK),   INTENT(IN)  :: lamb, t
      REAL(REALK),   INTENT(OUT) :: Fij(:,:)

      INTEGER(INTK), PARAMETER :: npoints = 10
      REAL(REALK)   :: rt(npoints), wt(npoints)
      REAL(REALK)   :: tmax, F(7)
      INTEGER(INTK) :: i,k, u,v, nt

      Fij = zero
      tmax = 6d0*lamb
      CALL fefmm_get_roots_and_weights(npoints,tmax,rt,wt)

      DO i = 1, FEdim 
         v = i + FEdim*((NFS-1)/2)  ! V-point

         DO k = 1, NFS*( (FEdim -1)/LIPN )   ! finite elements
            u = 1 + LIPN*(k-1)    ! basic rho-point

            F(:) = zero
            DO nt = 1, npoints
               ! 1st LIP
               F(1) = F(1) + wt(nt)*get_I1_LIP6(rt(nt),t,lamb,v,u)
               ! 2nd LIP
               F(2) = F(2) + wt(nt)*get_I2_LIP6(rt(nt),t,lamb,v,u)
               ! 3rd LIP
               F(3) = F(3) + wt(nt)*get_I3_LIP6(rt(nt),t,lamb,v,u)
               ! 4th LIP
               F(4) = F(4) + wt(nt)*get_I4_LIP6(rt(nt),t,lamb,v,u)
               ! 5th LIP
               F(5) = F(5) + wt(nt)*get_I5_LIP6(rt(nt),t,lamb,v,u)
               ! 6th LIP
               F(6) = F(6) + wt(nt)*get_I6_LIP6(rt(nt),t,lamb,v,u)
               ! 7th LIP
               F(7) = F(7) + wt(nt)*get_I7_LIP6(rt(nt),t,lamb,v,u)
            END DO

            DO nt = 0, 6
               Fij(u+nt,i) = Fij(u+nt,i) + F(nt+1)
            END DO

!if (abs(Fij(u,i)) > 1d-5) print *, "num_Fij: ", u,i, v-u, F(1)

         END DO

      END DO

   END SUBROUTINE fefmm_build_Fij_LIP6_numerical

!-------------------------------------------------------------------------------

   FUNCTION get_I1_LIP6(x,t,h,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: x, t
      REAL(REALK),   INTENT(IN) :: h
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_I1_LIP6

      INTEGER(INTK) :: k
      REAL(REALK)   :: f,g, fac

      k = v-u

      fac = one/(720d0*h**6)
      f = (x-h)* (x-2d0*h)* (x-3d0*h)* (x-4d0*h)* (x-5d0*h)* (x-6d0*h)
      g = t * (k*h - x)
      g = EXP(-g*g)
 
      get_I1_LIP6 = fac*f*g

   END FUNCTION get_I1_LIP6

!-------------------------------------------------------------------------------

   FUNCTION get_I2_LIP6(x,t,h,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: x, t
      REAL(REALK),   INTENT(IN) :: h
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_I2_LIP6

      INTEGER(INTK) :: k
      REAL(REALK)   :: f,g, fac

      k = v-u

      fac = -one/(120d0*h**6)
      f = x * (x-2d0*h)* (x-3d0*h)* (x-4d0*h)* (x-5d0*h)* (x-6d0*h)
      g = t * (k*h - x)
      g = EXP(-g*g)
 
      get_I2_LIP6 = fac*f*g

   END FUNCTION get_I2_LIP6

!-------------------------------------------------------------------------------

   FUNCTION get_I3_LIP6(x,t,h,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: x, t
      REAL(REALK),   INTENT(IN) :: h
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_I3_LIP6

      INTEGER(INTK) :: k
      REAL(REALK)   :: f,g, fac

      k = v-u

      fac = one/(48d0*h**6)
      f = x* (x-h)* (x-3d0*h)* (x-4d0*h)* (x-5d0*h)* (x-6d0*h)
      g = t * (k*h - x)
      g = EXP(-g*g)
 
      get_I3_LIP6 = fac*f*g

   END FUNCTION get_I3_LIP6

!-------------------------------------------------------------------------------

   FUNCTION get_I4_LIP6(x,t,h,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: x, t
      REAL(REALK),   INTENT(IN) :: h
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_I4_LIP6

      INTEGER(INTK) :: k
      REAL(REALK)   :: f,g, fac

      k = v-u

      fac = -one/(36d0*h**6)
      f = x* (x-h)* (x-2d0*h)* (x-4d0*h)* (x-5d0*h)* (x-6d0*h)
      g = t * (k*h - x)
      g = EXP(-g*g)
 
      get_I4_LIP6 = fac*f*g

   END FUNCTION get_I4_LIP6

!-------------------------------------------------------------------------------

   FUNCTION get_I5_LIP6(x,t,h,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: x, t
      REAL(REALK),   INTENT(IN) :: h
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_I5_LIP6

      INTEGER(INTK) :: k
      REAL(REALK)   :: f,g, fac

      k = v-u

      fac = one/(48d0*h**6)
      f = x* (x-h)* (x-2d0*h)* (x-3d0*h)* (x-5d0*h)* (x-6d0*h)
      g = t * (k*h - x)
      g = EXP(-g*g)
 
      get_I5_LIP6 = fac*f*g

   END FUNCTION get_I5_LIP6

!-------------------------------------------------------------------------------

   FUNCTION get_I6_LIP6(x,t,h,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: x, t
      REAL(REALK),   INTENT(IN) :: h
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_I6_LIP6

      INTEGER(INTK) :: k
      REAL(REALK)   :: f,g, fac

      k = v-u

      fac = -one/(120d0*h**6)
      f = x * (x-h)* (x-2d0*h)* (x-3d0*h)* (x-4d0*h)* (x-6d0*h)
      g = t * (k*h - x)
      g = EXP(-g*g)
 
      get_I6_LIP6 = fac*f*g

   END FUNCTION get_I6_LIP6

!-------------------------------------------------------------------------------

   FUNCTION get_I7_LIP6(x,t,h,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: x, t
      REAL(REALK),   INTENT(IN) :: h
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_I7_LIP6

      INTEGER(INTK) :: k
      REAL(REALK)   :: f,g, fac

      k = v-u

      fac = one/(720d0*h**6)
      f = x* (x-h)* (x-2d0*h)* (x-3d0*h)* (x-4d0*h)* (x-5d0*h)
      g = t * (k*h - x)
      g = EXP(-g*g)
 
      get_I7_LIP6 = fac*f*g

   END FUNCTION get_I7_LIP6

!-------------------------------------------------------------------------------

   FUNCTION get_F1_LIP6(t,lambda,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: lambda
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_F1_LIP6

      INTEGER(INTK) :: k
      REAL(REALK)   :: S1,S2, E1,E2, P,Q

      k = v-u

      E1 = exp(-t ** 2 * lambda ** 2 * (k - 0.6D1) ** 2)
      E2 = exp(-t ** 2 * k ** 2 * lambda ** 2)
      Q  = erf(t * k * lambda) - erf(-0.6D1 * lambda * t + t * k * lambda)

      S1 =  ((-k ** 5 / 0.1440D4 + 0.5D1 / 0.32D2 * k ** 2 - 0.137D3 / 0&
     &.720D3 * k + k ** 4 / 0.96D2 + 0.1D1 / 0.12D2 - 0.17D2 / 0.288D3 *&
     & k ** 3) / lambda / t ** 2 + (0.13D2 / 0.320D3 * k ** 2 - 0.317D3 &
     &/ 0.2880D4 * k + 0.11D2 / 0.120D3 - 0.7D1 / 0.1440D4 * k ** 3) / l&
     &ambda ** 3 / t ** 4 + (0.13D2 / 0.960D3 - 0.11D2 / 0.1920D4 * k) /&
     & lambda ** 5 / t ** 6)

      S2 =  ((-0.49D2 / 0.40D2 - 0.7D1 / 0.480D3 * k ** 4 + 0.35D2 / 0.2&
     &88D3 * k ** 3 - 0.49D2 / 0.96D2 * k ** 2 + k ** 5 / 0.1440D4 + 0.2&
     &03D3 / 0.180D3 * k) / lambda / t ** 2 + (0.175D3 / 0.576D3 * k + 0&
     &.7D1 / 0.1440D4 * k ** 3 - 0.21D2 / 0.320D3 * k ** 2 - 0.49D2 / 0.&
     &96D2) / lambda ** 3 / t ** 4 + (0.11D2 / 0.1920D4 * k - 0.7D1 / 0.&
     &240D3) / lambda ** 5 / t ** 6)

      P   = ((0.35D2 / 0.288D3 * k ** 4 - 0.49D2 / 0.40D2 * k + k ** 6 /&
     & 0.1440D4 + 0.203D3 / 0.180D3 * k ** 2 - 0.7D1 / 0.480D3 * k ** 5 &
     &- 0.49D2 / 0.96D2 * k ** 3 + 0.1D1 / 0.2D1) / t + (k ** 4 / 0.192D&
     &3 + 0.35D2 / 0.96D2 * k ** 2 - 0.7D1 / 0.96D2 * k ** 3 - 0.49D2 / &
     &0.64D2 * k + 0.203D3 / 0.360D3) / lambda ** 2 / t ** 3 + (-0.7D1 /&
     & 0.128D3 * k + 0.35D2 / 0.384D3 + k ** 2 / 0.128D3) / lambda ** 4 &
     &/ t ** 5 + 0.1D1 / lambda ** 6 / t ** 7 / 0.768D3)

      get_F1_LIP6 = S1*E1 + S2*E2 + ROOTPI*P*Q

   END FUNCTION get_F1_LIP6

!-------------------------------------------------------------------------------

   FUNCTION get_F2_LIP6(t,lambda,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: lambda
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_F2_LIP6

      INTEGER(INTK) :: k
      REAL(REALK)   :: S1,S2, E1,E2, P,Q

      k = v-u

      E1 = exp(-t ** 2 * lambda ** 2 * (k - 0.6D1) ** 2)
      E2 = exp(-t ** 2 * k ** 2 * lambda ** 2)
      Q  = erf(t * k * lambda) - erf(-0.6D1 * lambda * t + t * k * lambda)

      S1 =  ((k / 0.2D1 - 0.7D1 / 0.120D3 * k ** 4 - 0.77D2 / 0.120D3 *k&
     & ** 2 + k ** 5 / 0.240D3 + 0.71D2 / 0.240D3 * k ** 3) / lambda / t&
     & ** 2 + (-0.17D2 / 0.48D2 + 0.259D3 / 0.480D3 * k - 0.9D1 / 0.40D2&
     & * k ** 2 + 0.7D1 / 0.240D3 * k ** 3) / lambda ** 3 / t ** 4 + (-0&
     &.7D1 / 0.96D2 + 0.11D2 / 0.320D3 * k) / lambda ** 5 / t ** 6)

      S2 =  ((-0.31D2 / 0.48D2 * k ** 3 - k ** 5 / 0.240D3 + 0.3D1 + k *&
     &* 4 / 0.12D2 + 0.29D2 / 0.12D2 * k ** 2 - 0.87D2 / 0.20D2 * k) / l&
     &ambda / t ** 2 + (0.3D1 / 0.8D1 * k ** 2 - 0.155D3 / 0.96D2 * k - &
     &0.7D1 / 0.240D3 * k ** 3 + 0.29D2 / 0.12D2) / lambda ** 3 / t ** 4&
     & + (0.1D1 / 0.6D1 - 0.11D2 / 0.320D3 * k) / lambda ** 5 / t ** 6)

      P = ((k ** 5 / 0.12D2 - k ** 6 / 0.240D3 - 0.31D2 / 0.48D2 * k *&
     &* 4 + 0.29D2 / 0.12D2 * k ** 3 - 0.87D2 / 0.20D2 * k ** 2 + 0.3D1 &
     &* k) / t + (-0.31D2 / 0.16D2 * k ** 2 - k ** 4 / 0.32D2 + 0.5D1 / &
     &0.12D2 * k ** 3 - 0.87D2 / 0.40D2 + 0.29D2 / 0.8D1 * k) / lambda *&
     &* 2 / t ** 3 + (0.5D1 / 0.16D2 * k - 0.31D2 / 0.64D2 - 0.3D1 / 0.6&
     &4D2 * k ** 2) / lambda ** 4 / t ** 5 - 0.1D1 / t ** 7 / lambda ** &
     &6 / 0.128D3)

      get_F2_LIP6 = S1*E1 + S2*E2 + ROOTPI*P*Q

   END FUNCTION get_F2_LIP6

!-------------------------------------------------------------------------------

   FUNCTION get_F3_LIP6(t,lambda,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: lambda
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_F3_LIP6

      INTEGER(INTK) :: k
      REAL(REALK)   :: S1,S2, E1,E2, P,Q

      k = v-u

      E1 = exp(-t ** 2 * lambda ** 2 * (k - 0.6D1) ** 2)
      E2 = exp(-t ** 2 * k ** 2 * lambda ** 2)
      Q  = erf(t * k * lambda) - erf(-0.6D1 * lambda * t + t * k * lambda)

      S1 = ((-k ** 5 / 0.96D2 - 0.59D2 / 0.96D2 * k ** 3 + 0.13D2 / 0.9&
     &6D2 * k ** 4 + 0.107D3 / 0.96D2 * k ** 2 - 0.5D1 / 0.8D1 * k) / la&
     &mbda / t ** 2 + (0.33D2 / 0.64D2 * k ** 2 - 0.211D3 / 0.192D3 * k &
     &+ 0.7D1 / 0.12D2 - 0.7D1 / 0.96D2 * k ** 3) / lambda ** 3 / t ** 4&
     & + (0.31D2 / 0.192D3 - 0.11D2 / 0.128D3 * k) / lambda ** 5 / t ** &
     &6)

      S2 = ((-0.461D3 / 0.96D2 * k ** 2 + 0.137D3 / 0.96D2 * k ** 3 - 0&
     &.15D2 / 0.4D1 + 0.117D3 / 0.16D2 * k + k ** 5 / 0.96D2 - 0.19D2 / &
     &0.96D2 * k ** 4) / lambda / t ** 2 + (-0.461D3 / 0.96D2 + 0.685D3 &
     &/ 0.192D3 * k + 0.7D1 / 0.96D2 * k ** 3 - 0.57D2 / 0.64D2 * k ** 2&
     &) / lambda ** 3 / t ** 4 + (-0.19D2 / 0.48D2 + 0.11D2 / 0.128D3 * &
     &k) / lambda ** 5 / t ** 6)

      P = ((0.117D3 / 0.16D2 * k ** 2 - 0.461D3 / 0.96D2 * k ** 3 + k &
     &** 6 / 0.96D2 + 0.137D3 / 0.96D2 * k ** 4 - 0.19D2 / 0.96D2 * k **&
     & 5 - 0.15D2 / 0.4D1 * k) / t + (0.137D3 / 0.32D2 * k ** 2 + 0.117D&
     &3 / 0.32D2 + 0.5D1 / 0.64D2 * k ** 4 - 0.95D2 / 0.96D2 * k ** 3 - &
     &0.461D3 / 0.64D2 * k) / lambda ** 2 / t ** 3 + (0.137D3 / 0.128D3 &
     &- 0.95D2 / 0.128D3 * k + 0.15D2 / 0.128D3 * k ** 2) / lambda ** 4 &
     &/ t ** 5 + 0.5D1 / 0.256D3 / t ** 7 / lambda ** 6)

      get_F3_LIP6 = S1*E1 + S2*E2 + ROOTPI*P*Q

   END FUNCTION get_F3_LIP6

!-------------------------------------------------------------------------------

   FUNCTION get_F4_LIP6(t,lambda,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: lambda
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_F4_LIP6

      INTEGER(INTK) :: k
      REAL(REALK)   :: S1,S2, E1,E2, P,Q

      k = v-u

      E1 = exp(-t ** 2 * lambda ** 2 * (k - 0.6D1) ** 2)
      E2 = exp(-t ** 2 * k ** 2 * lambda ** 2)
      Q  = erf(t * k * lambda) - erf(-0.6D1 * lambda * t + t * k * lambda)

      S1 = ((0.49D2 / 0.72D2 * k ** 3 - k ** 4 / 0.6D1 - 0.13D2 / 0.12D2&
     & * k ** 2 + 0.5D1 / 0.9D1 * k + k ** 5 / 0.72D2) / lambda / t ** 2&
     & + (-0.13D2 / 0.24D2 + 0.173D3 / 0.144D3 * k + 0.7D1 / 0.72D2 * k &
     &** 3 - 0.5D1 / 0.8D1 * k ** 2) / lambda ** 3 / t ** 4 + (-0.3D1 / &
     &0.16D2 + 0.11D2 / 0.96D2 * k) / lambda ** 5 / t ** 6)

      S2 = ((-0.121D3 / 0.72D2 * k ** 3 + 0.31D2 / 0.6D1 * k ** 2 + 0.1&
     &0D2 / 0.3D1 - k ** 5 / 0.72D2 + k ** 4 / 0.4D1 - 0.127D3 / 0.18D2 &
     &* k) / lambda / t ** 2 + (0.31D2 / 0.6D1 + 0.9D1 / 0.8D1 * k ** 2 &
     &- 0.7D1 / 0.72D2 * k ** 3 - 0.605D3 / 0.144D3 * k) / lambda ** 3 /&
     & t ** 4 + (-0.11D2 / 0.96D2 * k + 0.1D1 / 0.2D1) / lambda ** 5 / t&
     & ** 6)

      P = ((-0.121D3 / 0.72D2 * k ** 4 + 0.10D2 / 0.3D1 * k - k ** 6 /&
     & 0.72D2 + k ** 5 / 0.4D1 + 0.31D2 / 0.6D1 * k ** 3 - 0.127D3 / 0.1&
     &8D2 * k ** 2) / t + (0.5D1 / 0.4D1 * k ** 3 - 0.5D1 / 0.48D2 * k *&
     &* 4 + 0.31D2 / 0.4D1 * k - 0.121D3 / 0.24D2 * k ** 2 - 0.127D3 / 0&
     &.36D2) / lambda ** 2 / t ** 3 + (-0.5D1 / 0.32D2 * k ** 2 + 0.15D2&
     & / 0.16D2 * k - 0.121D3 / 0.96D2) / lambda ** 4 / t ** 5 - 0.5D1 /&
     & 0.192D3 / t ** 7 / lambda ** 6)

      get_F4_LIP6 = S1*E1 + S2*E2 + ROOTPI*P*Q

   END FUNCTION get_F4_LIP6

!-------------------------------------------------------------------------------

   FUNCTION get_F5_LIP6(t,lambda,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: lambda
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_F5_LIP6

      INTEGER(INTK) :: k
      REAL(REALK)   :: S1,S2, E1,E2, P,Q

      k = v-u

      E1 = exp(-t ** 2 * lambda ** 2 * (k - 0.6D1) ** 2)
      E2 = exp(-t ** 2 * k ** 2 * lambda ** 2)
      Q  = erf(t * k * lambda) - erf(-0.6D1 * lambda * t + t * k * lambda)

      S1 = ((-k ** 5 / 0.96D2 - 0.5D1 / 0.16D2 * k + 0.11D2 / 0.96D2 * k&
     & ** 4 + 0.61D2 / 0.96D2 * k ** 2 - 0.41D2 / 0.96D2 * k ** 3) / lam&
     &bda / t ** 2 + (-0.145D3 / 0.192D3 * k + 0.7D1 / 0.24D2 - 0.7D1 / &
     &0.96D2 * k ** 3 + 0.27D2 / 0.64D2 * k ** 2) / lambda ** 3 / t ** 4&
     & + (0.23D2 / 0.192D3 - 0.11D2 / 0.128D3 * k) / lambda ** 5 / t ** &
     &6)

      S2 = ((0.107D3 / 0.96D2 * k ** 3 - 0.307D3 / 0.96D2 * k ** 2 - 0.&
     &15D2 / 0.8D1 + k ** 5 / 0.96D2 + 0.33D2 / 0.8D1 * k - 0.17D2 / 0.9&
     &6D2 * k ** 4) / lambda / t ** 2 + (0.7D1 / 0.96D2 * k ** 3 - 0.307&
     &D3 / 0.96D2 - 0.51D2 / 0.64D2 * k ** 2 + 0.535D3 / 0.192D3 * k) / &
     &lambda ** 3 / t ** 4 + (-0.17D2 / 0.48D2 + 0.11D2 / 0.128D3 * k) /&
     & lambda ** 5 / t ** 6)

      P = ((-0.17D2 / 0.96D2 * k ** 5 - 0.15D2 / 0.8D1 * k - 0.307D3 /&
     & 0.96D2 * k ** 3 + k ** 6 / 0.96D2 + 0.33D2 / 0.8D1 * k ** 2 + 0.1&
     &07D3 / 0.96D2 * k ** 4) / t + (0.33D2 / 0.16D2 - 0.307D3 / 0.64D2 &
     &* k + 0.107D3 / 0.32D2 * k ** 2 - 0.85D2 / 0.96D2 * k ** 3 + 0.5D1&
     & / 0.64D2 * k ** 4) / lambda ** 2 / t ** 3 + (0.107D3 / 0.128D3 - &
     &0.85D2 / 0.128D3 * k + 0.15D2 / 0.128D3 * k ** 2) / lambda ** 4 / &
     &t ** 5 + 0.5D1 / 0.256D3 / t ** 7 / lambda ** 6) 

      get_F5_LIP6 = S1*E1 + S2*E2 + ROOTPI*P*Q

   END FUNCTION get_F5_LIP6

!-------------------------------------------------------------------------------

   FUNCTION get_F6_LIP6(t,lambda,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: lambda
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_F6_LIP6

      INTEGER(INTK) :: k
      REAL(REALK)   :: S1,S2, E1,E2, P,Q

      k = v-u

      E1 = exp(-t ** 2 * lambda ** 2 * (k - 0.6D1) ** 2)
      E2 = exp(-t ** 2 * k ** 2 * lambda ** 2)
      Q  = erf(t * k * lambda) - erf(-0.6D1 * lambda * t + t * k * lambda)

      S1 = ((-k ** 4 / 0.24D2 + 0.7D1 / 0.48D2 * k ** 3 - 0.5D1 / 0.24D2&
     & * k ** 2 + k ** 5 / 0.240D3 + k / 0.10D2) / lambda / t ** 2 + (0.&
     &127D3 / 0.480D3 * k - 0.3D1 / 0.20D2 * k ** 2 + 0.7D1 / 0.240D3 * &
     &k ** 3 - 0.17D2 / 0.240D3) / lambda ** 3 / t ** 4 + (-0.19D2 / 0.4&
     &80D3 + 0.11D2 / 0.320D3 * k) / lambda ** 5 / t ** 6)

      S2 = ((-0.27D2 / 0.20D2 * k + k ** 4 / 0.15D2 + 0.3D1 / 0.5D1 - 0&
     &.19D2 / 0.48D2 * k ** 3 + 0.13D2 / 0.12D2 * k ** 2 - k ** 5 / 0.24&
     &0D3) / lambda / t ** 2 + (-0.95D2 / 0.96D2 * k + 0.3D1 / 0.10D2 * &
     &k ** 2 + 0.13D2 / 0.12D2 - 0.7D1 / 0.240D3 * k ** 3) / lambda ** 3&
     & / t ** 4 + (0.2D1 / 0.15D2 - 0.11D2 / 0.320D3 * k) / lambda ** 5 &
     &/ t ** 6)

      P = ((0.13D2 / 0.12D2 * k ** 3 + k ** 5 / 0.15D2 - 0.19D2 / 0.48&
     &D2 * k ** 4 - k ** 6 / 0.240D3 - 0.27D2 / 0.20D2 * k ** 2 + 0.3D1 &
     &/ 0.5D1 * k) / t + (-k ** 4 / 0.32D2 - 0.19D2 / 0.16D2 * k ** 2 - &
     &0.27D2 / 0.40D2 + 0.13D2 / 0.8D1 * k + k ** 3 / 0.3D1) / lambda **&
     & 2 / t ** 3 + (-0.19D2 / 0.64D2 - 0.3D1 / 0.64D2 * k ** 2 + k / 0.&
     &4D1) / lambda ** 4 / t ** 5 - 0.1D1 / t ** 7 / lambda ** 6 / 0.128&
     &D3)

      get_F6_LIP6 = S1*E1 + S2*E2 + ROOTPI*P*Q

   END FUNCTION get_F6_LIP6

!-------------------------------------------------------------------------------

   FUNCTION get_F7_LIP6(t,lambda,v,u)

      IMPLICIT NONE
      REAL(REALK),   INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: lambda
      INTEGER(INTK), INTENT(IN) :: v  ! this is the V-point
      INTEGER(INTK), INTENT(IN) :: u  ! this is the rho-point
      REAL(REALK) :: get_F7_LIP6

      INTEGER(INTK) :: k
      REAL(REALK)   :: S1,S2, E1,E2, P,Q

      k = v-u

      E1 = exp(-t ** 2 * lambda ** 2 * (k - 0.6D1) ** 2)
      E2 = exp(-t ** 2 * k ** 2 * lambda ** 2)
      Q  = erf(t * k * lambda) - erf(-0.6D1 * lambda * t + t * k * lambda)

      S1 = ((-k ** 5 / 0.1440D4 - 0.1D1 / 0.12D2 + k ** 4 / 0.160D3 - k &
     &/ 0.36D2 + 0.13D2 / 0.480D3 * k ** 2 - 0.31D2 / 0.1440D4 * k ** 3)&
     & / lambda / t ** 2 + (-0.119D3 / 0.2880D4 * k - 0.7D1 / 0.1440D4 *&
     & k ** 3 + 0.7D1 / 0.320D3 * k ** 2) / lambda ** 3 / t ** 4 + (0.1D&
     &1 / 0.192D3 - 0.11D2 / 0.1920D4 * k) / lambda ** 5 / t ** 6)

      S2 = ((0.137D3 / 0.720D3 * k - k ** 4 / 0.96D2 + 0.17D2 / 0.288D3&
     & * k ** 3 - 0.5D1 / 0.32D2 * k ** 2 - 0.1D1 / 0.12D2 + k ** 5 / 0.&
     &1440D4) / lambda / t ** 2 + (0.7D1 / 0.1440D4 * k ** 3 - 0.5D1 / 0&
     &.32D2 + 0.85D2 / 0.576D3 * k - 0.3D1 / 0.64D2 * k ** 2) / lambda *&
     &* 3 / t ** 4 + (0.11D2 / 0.1920D4 * k - 0.1D1 / 0.48D2) / lambda *&
     &* 5 / t ** 6)

      P = ((-0.5D1 / 0.32D2 * k ** 3 - k ** 5 / 0.96D2 + 0.17D2 / 0.28&
     &8D3 * k ** 4 + 0.137D3 / 0.720D3 * k ** 2 - k / 0.12D2 + k ** 6 / &
     &0.1440D4) / t + (-0.5D1 / 0.96D2 * k ** 3 - 0.15D2 / 0.64D2 * k + &
     &0.17D2 / 0.96D2 * k ** 2 + k ** 4 / 0.192D3 + 0.137D3 / 0.1440D4) &
     &/ lambda ** 2 / t ** 3 + (k ** 2 / 0.128D3 - 0.5D1 / 0.128D3 * k +&
     & 0.17D2 / 0.384D3) / lambda ** 4 / t ** 5 + 0.1D1 / t ** 7 / lambd&
     &a ** 6 / 0.768D3)

      get_F7_LIP6 = S1*E1 + S2*E2 + ROOTPI*P*Q

   END FUNCTION get_F7_LIP6

!-------------------------------------------------------------------------------

END MODULE fefmm_LIP6_Fij
